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Abstract 

The multiscaling properties of the mixed Obukhov-Novikov shell model of 
turbulence are investigated numerically and compared with those of the com- 
plex GOY model, mostly studied in the recent years. Two types of generic 
singular fluctuations are identified : first, self-similar solutions propagating 
from large to small scales and building up intermittency, second, complex 
time singularities inhibiting the cascade and promoting chaos. A simple and 
robust method is proposed to track these objects. It is shown that the scaling 
exponent of self-similar solutions selected by the dynamics is compatible with 
large order statistics whenever it departs enough from the Kolmogorov value. 
Complex time singularities on the other hand get trapped on the last shells, 
when the proportion of Novikov interactions exceeds a critical value which is 
argued to mark the boundary between chaotic and regular dynamics in the 

limit of infinite Reynolds number. 
PACS numbers 03.40.G, 47.27 



Typeset using REVTgX 



1 



I. INTRODUCTION 



The shell models of turbulence have recently attracted a lot of interest as a useful tool 
for mimicking the Navier-Stokes dynamics. In the simplest scalar models, one places veloc- 
ity variables on a one-dimensional array of wavevectors of the form k n = koQ n , where the 
integer n labels the shell from to, ideally, +00 and Q is a scale parameter fixing the step of 
the cascade. The time evolution of shell-velocities is governed by ordinary differential equa- 
tions with quadratic non-linearities, whose strength grows like k n , deterministic forcing at 
large scales and viscous dissipation at small ones. The couplings between shells are usually 
local and chosen in such a way that the total kinetic energy is conserved in the absence of 
forcing and viscous effects. These hydrodynamic systems display strong departure from the 
naive scaling expected on the basis of Kolmogorov- like dimensional analysis, which shows 
up in particular in the higher order moments of velocity. Following the seminal work of 
Okhitani and Yamada |l|,@], Jensen, Paladin and Vulpiani || found in a particular shell 
model, nowadays referred to as the GOY model, multiscaling properties very close to those 
of real turbulent flows. Most subsequent studies in this field have therefore concentrated on 
the GOY model and important progress was made towards a deeper understanding of its 
behaviour in recent publications (we shall go back to some of these results in the bulk 
of the paper). 

The GOY model uses complex velocity variables and interactions among all triads made 
up of three different neighbouring shells. One may wonder whether these two features are 
necessary to produce "good" chaotic properties. In order to clarify this question, we report 
in this paper a mostly numerical investigation of the scaling properties of a simpler class 
of shell models, which results from the linear superposition of two chains introduced in the 
early 70's by Obukhov and Novikov || (for an historical insight into the field and a 
comprehensive review of the huge russian litterature concerned with cascade-like systems 
under various disguises, see for instance 0). In the Obukhov- Novikov model, hereafter 
referred to as the ON model, the velocity variables are real and interacting triads involve 
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only two neighbouring shells. The structure of non-linearities depends as in the GOY model 
on a single parameter (together with the scale parameter Q), which fixes in that case the 
relative proportion of the two basic chains. Both models exhibit qualitatively the same 
phenomenology. When the proportion of Novikov-like interactions (favouring the transfer 
of energy towards small scales) is high enough, the system relaxes to a time independent 
state with Kolmogorov scaling properties. As Obukhov-like interactions (favouring on the 
contrary the backflow of energy towards large scales) take over, the system transits through 
a Ruelle-Takens scenario into a chaotic state with stochastic fluctuations. There is clearly 
multifractality close to the transition, even if it looks less pronounced than the one observed 
in the GOY model for usual values of parameters. 

We switch to more deterministic concerns in the second part of this paper, which aims 
at characterizing singular fluctuations able to form in the ON model or more generally in 
any one-dimensional shell model. We shall argue that self-similar or soliton-like solutions 
of the equations of motion in the inertial range are the building blocks of intermittency, 
while movable singularities occuring at complex times induce chaos by inhibiting the energy 
cascade. We are not aware of any previous study of the structure of complex time singu- 
larities in shell models. In contrast, self-similar solutions were already considered within 



the context of the ON model by Siggia |10[ and later on in more details by Nakano |TT 
They are curiously absent of more recent works. We propose here an efficient method for 
identifying such solutions without any a priori assumption on their shape. We find that 
the set of dynamically accessible self-similar solutions is in fact limited to one single object 
(as Nakano's results suggested it). This proves that multiscaling properties should not be 
ascribed to the existence of a large manifold of singular behaviours. The exponent z con- 
trolling the multiplicative growth of these particular solutions, accounts in a satisfying way 
for the asymptotic scaling properties of high order velocity or energy transfer moments. It 
is easy to extend this analysis to the GOY model, where basically the same conclusions 
concerning the unicity of solutions can be drawn. However self-similar solutions in that case 
are very mild and do not seem to play a major role in the statistics at high orders. 



The paper is organized as follows : in Section 2 we specify the conventions used in our com- 
putations for normalizing variables and parameters and describe some general properties of 
the ON model. Section 3 presents statistical results obtained from numerics. The emphasis 
is put on scaling exponents of the moments of energy transfer and their evolution with the 
relative proportion of Obukhov and Novikov interactions. Our goal here is to provide the 
reader with data and facts, disentangled from any theoretical interpretation. An attempt of 
comparison with the GOY model is made. Section 4 is devoted to the hunt for self-similar 
solutions and a confrontation of their scaling properties with the statistics of the model 
at large orders. Complex time singularities are introduced and studied in Section 5, while 
perspectives and conclusions are briefly outlined in Section 6. 



II. GENERAL PROPERTIES OF THE MODEL 

A. Definitions and basic considerations 

As already said in the Introduction, scalar shell models define a velocity variable u n , real 
or complex, on a one-dimensional array of wave-vectors k n = k Q n where the integer n runs 
from to N. In most of the paper we shall restrict ourselves to the case of real variables. 
It simplifies notations to consider that the u n form a (N + l)-dimensional vector u. The 
equation of motion then takes the following form 

^u = N[u]+F-D (1) 

where the three vectors N, F, D embody respectively the non-linearities, the external forcing 
and the dissipation. We only considered a deterministic forcing acting on the zero*' 1 shell 
and usual viscous dissipation, which means 

F n = f5 n ,o , D n = vk 2 n u n (2) 

where v is the kinematic viscosity. The nonlinear kernel N is quadratic in the u n , with a 
coupling constant growing like k n in order to reproduce the hierarchy of characteristic times 
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1 N 

of the Navier-Stokes dynamics. It must also conserve the total kinetic energy E = — ^ u n . 

^ n=0 

If interactions between shells (which are always supposed to be local) do not extend beyond 
nearest neighbours, the most general expression for the n th component of iV is 

N n [u] = aQ^ [knUn-tUn - k n+1 u 2 n+1 ] + fi\k n u 2 n _ x - k n+l u n u n+l ] (3) 

(this formula remains valid on the two boundaries n = and n = N, provided it_i = u^ + \ = 
is assumed). 

The model appears like the linear superposition of the Obukhov-Gledzer (OG) and Novikov- 
Desniansky (ND) chains, with respective weig hts aQ 2/3 and (3. We shall assume a, (3 > 
and, without loss of generality, a + f3 — 1. Since, on the average, the u n decrease like fc~ 1//3 
according to Kolmogorov-scaling, it is convenient to introduce a new set of variables 4> n by 
the relation 

U n = Q~^<Pn (4) 

The equations for read 

j t $=k Q-3N[$]+F-D (5) 
where the expression of the n th component of iV is now 

N n [$\ = Q 2 ? [(a0 n _x0 n + /ty^) - (a4> 2 n+1 + /30 n n+1 )] (6) 

We still have the freedom to set to unity the forcing amplitude and the coefficient in front 
of N by non-dimensionalizing in the proper way time and velocities. The final form of the 
equations (as they were used in the numerical investigations reported in Section 3) is 

^ n = N n [^]+6 nfi -^Q 2n ^ n (7) 
at ri 



The Reynolds number R has been defined as R = — > 

v 



\ 3 . Equations @ and (0) make 



energy conservation quite obvious. Indeed, in the limit R = +oo and for n > 1, the energy 
E n = — = -^<P n Q~~ carried by the n th shell obeys the equation 



~^^ n — 6 n — e n+l (8) 

where 

e n = <fi n -l<fin{a<fin + P<t>n-l) ( 9 ) 

is the energy flux from the (n — l) th to the n th shell. Kolmogorov scaling corresponds to 
(fi n = C te or more fundamentally to e n = C te throughout the cascade. 

The physics of the model, as defined by equation (0), depends on three parameters, 
namely : the step of the cascade Q, the proportion of Novikov interactions /3, and the 
Reynolds number R. The number of shells will not matter, provided the truncation is done 
far beyond the Kolmogorov dissipative scale, where viscous effects become of the same order 

as inertial ones. Assuming <fi n = 0(1), the index of the dissipative shell is given by the 

Q 2Nd ™& 

condition — — ~ Q » . One should however pay attention to the fact that the stronger the 
R 

fluctuation, the smaller the scale at which it will be effectively dissipated. Since <fi n can grow 

at most like Q~s (this corresponds to the extreme case of a fluctuation carrying a constant 

4 

energy through the cascade), we conclude that iV should be an integer between rid and -n^, 
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where we have defined as 

3Log# 
4LogQ 

In our numerical study of the statistical properties of the model, we took Q = 2, R = 10 5 
and let vary f3 between and 1. The choice of 18 shells (N =17) turned out to ensure the 
absence of any spurious boundary effect. 

Since we shall allude sometimes to the complex GOY model, we close this section by 
writing down the version of it we used in our computations. With complex variables <fi n 
rescaled in the way described just before, the equations read : 

± ( fi = (2-e)[l-m]-^ (Ha) 



^i = Q*0S&-0S#]-^i (lib) 
-<f> n = [(1 - e)ft_ift_a + e<f>U<f>* n+ i ~ ft+iftfJ ~ \K 2<n<N (11c) 



>N+\ — <PN+2 



(lid) 



Here again the structure of non-linearities is fixed by a single parameter e, varying as (3 
between and 1. We put for simplicity the forcing on the zero th shell, as first tried in 
||. The factor —i, which is usually kept in front of the non-linear terms, as a remnant 
of the Navier-Stokes equation, has been absorbed in an innocuous redefinition of variables 
{u n — > iu n ). The forcing also can be assumed to be real (here F n = (2 — e)5 n fi), without loss 
of generality, thanks to the invariance of equations under the following phase transformation 
(see for instance H) 

03n -> e^03n (12a) 
3 n+l -> e iei 03n+l (12b) 
03n+2 -> e- i(e ° +9l V3n + 2 (12c) 

where an d ^1 are arbitrary angles. There are obvious differences between the real ON 
model and the complex GOY model, which make their comparison interesting in its own 
right. The phase space which may be explored by the GOY model is a priori larger, at 
least for initial conditions not purely real. Also the range of non-linear interactions is wider 
than in the first model. Lastly, the GOY model admits a second quadratic invariant besides 
kinetic energy, which is thought to play a prominent role in fixing its statistical properties 



12]. Such an extra invariant is definitely absent in the ON model. 



B. Fixed points and qualitative description of the phase diagram 



First of all, let us say a few words about the existence and the nature of the fixed points 
of the model, which lead to Kolmogorov scaling. In the absence of dissipation and for an 
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infinite number of shells, <\> n = 1 is an obvious solution for any value of (3. However, once 
viscous dissipation is introduced, it is easy to see that static solutions making physical sense 
can exist only for (3 ^ 0. Indeed for a = 1 (OG chain), one has to solve for every n > 

Q^(0 n -10»-#+l) = ^n (13) 

In order to balance the leading order terms of this equation for n — > +oo, n must behave 

Q 3 

like — — which is not acceptable. This absence of any fixed point in the presence of 
dissipation for the Obukhov chain is consistent with the dynamical behaviour one observes 
in this case. A solitary wave of negative amplitude appears invariably at the end of the 
cascade, which carries back energy towards large scales. In the final state, energy oscillates 
back and forth between the first two shells, the other ones being inactive. 



Things are different, as soon as a finite amount of Novikov interactions is introduced. 
This is because the term (fi^-i i n equation (|^) favours energy transfer to small scales. At a 
static level, assuming a rapid decay of the spectrum on the ultraviolet side (i.e. n _i >> 
<p n » <f> n +i), one has now to achieve the balance 

0±2 Q 3 - 



R 



(14) 



4{n + 2) 

This yields the following general solution <p n = — — — exp(— b2 n ) with b > otherwise 

pR 

arbitrary. More physically, we may rewrite <p n in the dissipative range as 

4(n-n d +2) 

(p n = ^— j exp(-a2 (n -" d) ) for n > n d (15) 

with rid given by fllPl). These considerations suggest the following procedure for computing 
the fixed point § e n in the presence of dissipation. The condition of equilibrium for the n th 
shell being quadratic in n _i, it can be used to express <t> n -\ in terms of 4> n and 4> n+ ±. Going 
from the last shell with <pN parameterized as in fllBD to the first one, and keeping at each 
step the positive root, one ends up with N + 1 numbers 0jv(a), . . . , 4>i(a), 4>o(a). The last 
equation on the zero t/l shell yields then a solvability condition for the parameter a 

8 



F(a) = 1 - (a$(a) + p(p (a)M a )) = (16) 

Figures [I] and |2] show how things work for Q = 2, R = 10 5 and (3 = 0.348. At least 
in the range of Reynolds numbers investigated here, there exists only one solution and the 
resulting fixed point is drawn on Fig|| The emergence of oscillations of period 2 and growing 
amplitude on the infrared side of the spectrum, is a generic feature of the ON model. Similar 
phenomena involving period 3 are known to plague the GOY model. Here the explanation 
is the following : in the inertial range (i.e. 1 < n <C rid), the equation to be solved reduces 
to 

q-\l + Cq- 1 ) - q n+l {q n+l + C) = (17) 

where one has introduced the ratio q n = -^-^ and C = Equation ( |T7D defines a map 
Qn — 9(<ln+i), whose fixed point q n = 1 is easily seen to be unstable for C < 1. Only for 
discrete values of the Reynolds number (such that one of the crossing points in Fig.|I] has 
ordinate 1), can one hope to get rid of these oscillations. We should say that this odd-even 
disymmetry, though a pathology of the ON model, is less visible in the chaotic state to be 
described below. Furthermore, it does not affect other physical quantities like the energy 
transfer e n defined in (|9|). 

We turn now to a qualitative discussion of the phase diagram of the model, which is 
observed as (3 varies (for the particular values Q = 2 and R = 10 5 ). Although the static 
solution exists for every strictly positive /3, it becomes unstable for (3 < 0.355 ± 10~ 3 and 
evolves towards a periodic limit cycle through a first Hopf bifurcation. A scenario a la Ruelle 
Takens, similar to the one already discovered in the GOY model 0, leads then to chaos for 
(3 < (3* = 0.349 ± 10 -3 . We did not try to get very precise estimates of these two thresholds, 
whose position is expected to vary with R. From a slightly unrelated analysis of movable 
singularities of the ON model, presented in Section 5, we shall speculate later in the paper 
that the asymptotic value of (3* in the limit of infinite Reynolds number is of order 0.394. 
In the quasiperiodic regime, the shells oscillate in a coherent way around the static fixed 
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point discussed above. Their oscillations remain of moderate amplitude even close to the 
transition, and one has with a very good accuracy, for all integers n and p, (<f>^(t)) = (<^) p 
(where (. . .) denotes a temporal average). 

The behaviour in the chaotic phase is rather simple to understand, far from the transition. 
For values of (3 not greater than say 0.28, the dynamics consists in well isolated pulses emit- 
ted from the forced shell, after it has reached a significant level. The pulse propagates down 
the scales, through almost inactive shells, leaving behind a finite amount of energy. After 
being stopped by dissipation, it gives rise to a rather well characterized pulse of negative 
amplitude, which carries back to the large scales most of the energy that dissipation failed 
to absorb. This is quite close to what is observed in the Obukhov limit, although here things 
repeat themselves in a slightly disordered manner. The picture gets more complicated as (3 
gets closer to (3* (typically for 0.30 < (3 < 0.349). Shells in the inertial range remain almost 
always active and they form a noisy background out of which only the biggest fluctuations 
develop like singular pulses. Negative excursions of variables, triggering backflow of energy, 
occur very rarely. Splitting of pulses as well as partial reflection of some of them may be 
observed all the way along the cascade. The next Section will help to quantify a bit these 
statements. 



III. MULTISCALING PROPERTIES OF THE CHAOTIC PHASE 

Moments of the variables 0„, defined as S n ^ p = (4%), are good tools for characterizing 
the intermittency of shell models. According to (f|), exponential growth of any of the S^p's 
with n in the inertial range is the sign of deviations from K41-scaling. However we shall 
rather concentrate on moments of the energy flux e n defined in ([|), namely on quantities 
S n .p = (l^nl 3 ) The reason for this choice has already been discussed in || within the context 
of the GOY model. The point is that for moderate p, the oscillations we found in the static 
solution still contaminate S ntP . Although they are much smaller, they prevent us from an 
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accurate determination of the scaling exponents. As far as scaling properties are concerned, 
the E n p 's provide a valuable alternative, because they are free from any spurious oscillations. 
On the other hand one must be aware that the absolute value introduces a considerable bias 
for the lowest moments if the energy cascade has a poor positive "smooth component", as it 
is the case far away from the transition where the mean energy transfer is very small. When 
on the contrary the smooth component is important, e n is almost always positive and (|e„|) 
is very close to a constant in the inertial range as (e n ) should be in a statistically stationnary 
state. In any case the bias disappears at higher values of p since the strongest fluctuations 
are always positive. 

We have calculated S n p for p up to 12. For each run we integrated the equations over 1000 
turn-over or unit times and sampled the signal with a step much smaller than the character- 
istic time-scales of the more intense structures (in practice we took At = 10~ 4 ). The length 
of each run was enough to ensure the stationarity of the statistics, at least for moments of 
order p smaller than 7 — 8. As shown on Fig. ^ for (3 = 0.33, these quantities obey nice 
power laws on a rather wide range of shell numbers. We may thus define exponents a p such 
that in the inertial range E HiP ~ k n ap (the relation between a p and the usual exponent ( p 
associated to the velocity field is : ( p = | — a p ). 

To extract values of the scaling exponents a p from our data, we again followed the pro- 
cedure outlined in (though not with the same refinement!). The inertial or fitting range 
was determined as the interval of values of n, for which a least square fit of the data to a 
straight line, give o\ the closest to zero. A fitting range 5 < n < 10 proved to be the best 
for all values of (3 we looked at, with o\ as small as 0.01 for f3 ~ 0.33. Actually the highest 
moments (p > 7) allow a wider range for an estimation of the exponents, whereas the lowest 
depend only slightly on this length. In order to get a rough estimate of statistical errors 
due to the finite length of our temporal signal, we repeated the same operations many times 
(typically 5), taking as new initial conditions the state obtained at the end of the precedent 
run. Error bars drawn in the next figures take into account this only source of uncertainty. 
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Our results for the scaling exponents a p are summarized in Fig. |5| for three typical values of 
(3 : (3 = 0.28, = 0.33, and (3 = 0.343. 

As concerning the scattering of data, we see that it remains small for p < 8, except when we 
get very close to the transition, a not too surprising fact. The overall shape of the <r p -curve 
illustrates the distinction made at the end of Section 2 between two chaotic regimes. For 
(3 = 0.28, i.e. rather far from the transition, one crosses over rapidly (p > 5) towards a 
linear growth of exponents with p. The exponents take rather big negative values for p < 3, 
presumably because of the important role played by energy backflows in this case. For the 
two other values of (3 we investigated, the cross-over region is significantly wider. There is 
clearly curvature and henceforth multifractality, even if for p > 7 — 8 an asymptotic linear 
regime ultimately sets in. Figure |6| offers a magnified view on the cross-over region. The 
fact that exponents still vary strongly between (3 = 0.33 and (3 = 0.343 proves the influence 
of the proximity of the transition on the physics probed by the moments of corresponding 
order. On the other hand, it can be noted that the asymptotic slope remains almost the 
same. This feature will be explained in the next Section in terms of self-similar solutions 
parameterizing the biggest fluctuations of our hydrodynamic system. 

As a way of checking our numerics, and also in order to get some insight into the differ- 
ences between the two models, we have performed strictly analogous computations on the 
complex GOY model for two different values of e (with Q = 2, R = 10 5 and iV = 17). The 
first one, e = 0.75, lies far from the transition, which was found in || to occur at e* = 0.395.. 
for a Reynolds number roughly the halfth of ours. The second one is the standard e = 0.5, 
known to lead to scaling properties in good agreement with experiment || . Figure [7| shows 
on the same graph the scaling exponents for the ON model with (3 = 0.28 and the GOY 
one for e = 0.75. The resemblance is striking, apart from from slightly different asymptotic 
slopes of the two curves and a cross-over region a bit wider for the GOY model. In Fig. |8|, 
the same comparison is made between the two models, this time for (3 = 0.33 and e = 0.5. 
Now differences show up, especially at large orders. While a linear monofractal behaviour 
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has definitely set in for p > 7 in the case of the ON model, the local slope of the a p does 
not stop increasing in the GOY model, even though this is hardly perceptible to the eye. 
Also significant is the lack of convergence for p > 8. We note that the points marking the 
upper error bars, which were obtained from one particular run among six of equal temporal 
length, are surprisingly well fitted by the formula proposed by She and Leveque |TB|] for real 
Navier-Stokes turbulence 

*P = !-2[l-(f)«] (18) 

up to the highest values of p investigated here. A physical reason for this good agreement 
observed by many others is missing. We did not try to tighten our error bars by increasing 
the length of numerical integration since this was beyond the scope of this work. The results 
of the next Section will somehow confirm the peculiar nature of intermittency in the GOY 
model for such values of parameters as e = 0.5 and Q = 2. 



IV. STUDY OF SELF-SIMILAR SOLUTIONS 

Let us take for granted from the results of the preceding Section that the scaling expo- 
nents <j p grow asymptotically like 7p at large p, with 7 a positive number depending on the 
parameter (3 only. This means that the amplitude of fluctuations carrying the system away 
from the K41 fixed point cannot grow from shell to shell more rapidly than Q in . On the 
other hand, the fact that 7 takes a finite value even close to the transition suggest that such 
fluctuations are efficient as soon as the instability threshold is passed. It is the purpose of 
this Section to identify the set of singular fluctuations that the ON model can admit. 
Since we are now interested in nonlinear instabilities occuring in the inertial range, we may 
forget about forcing and dissipation, and think of the shell number n as running from —00 

In 

to +00. Let us rewrite the equation of motion in terms of new variables b n = Q n u n = Q~(f> n 
{b n is nothing dimensionally but the gradient of the velocity field). We get from (g) and (0), 

2 

after absorbing the factor aQa into a rescaling of time : 
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jb n = N n [b] = (6 n _ x 6 n + £g« 6^0 - ^(ft^i + ^QS&A+i) (19) 



Since iV[fc] does not depend explicitely on n and is quadratic in b, the set of equations (1191), 
for — oo < n < +00, support formally self-similar solutions of the type : 

= ^ /(g " 2(r " t)] = Q nz 3(Q nz ^ - *)) (20) 

In the equation above, t* is the critical time at which, in the absence of dissipation, the 
fluctuation reaches the end of the cascade. The scaling exponent z is a priori arbitrary. 
However z = | gives back Kolmogorov scaling, while z = 1 corresponds to the extreme 
situation of a fluctuation carrying a constant energy. One expects therefore | < z < 1 on 
physical grounds. 

Self-similar solutions, if they exist, are obviously good candidates for describing the growth 
of singular fluctuations. The question then arises whether many values of z are dynamically 
accessible, which would be a natural source of multifractality, or whether on the contrary 
a single z is selected. In that case, one should check that 7 = (z — |), since moments S njP 
are dominated by extreme fluctuations for high values of p. Self-similar solutions together 



with their exponent z have already been determined by Nakano for the ON model ||11| . He 
used a rather cumbersome iterative method to find them and we were not convinced he 
had exhausted the whole set of possibilities in his work. This is why we came back to this 
problem and were led to develop a procedure to be described below, which is quite efficient 
and easily extended to any shell model. It should however be said from the beginning that 
our results about the ON model are in complete agreement with the conclusions reached in 

0- 



By plugging the Ansatz ( p0|) into ( ]19"1) and introducing the logarithmic variable £ = 
n H — p— ^ log(t* — t), one arrives at the following equation for / (/ is actually divided by 
z log Q to make the result a bit simpler) 
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no -ziogQfio = (/(^ - i)/(o + - 1) 2 ) - i(/ce + 1) 2 + -Qhiofn + 1)) 

(21) 

If square integrability of / is required, Eq. (|2lD is nothing but a non-linear eigenvalue problem 
for the unknown z, which is very difficult to solve directly, either analytically or numerically. 
To make progress, we can try to approach / dynamically. Rather than coming back to the 
original equations of the model, let us introduce a fictitious dynamics leaving the norm of 
the (N + l)-dimensional vector b invariant 

d 7 -r-, <N\b],b> r , . 

—b = N[b\ b 22 

dr <b,b> 

In the equation above, N[b] is the vector of components iV n [&], whose expression was given 
in Eq. flI5|) and (A,B) = Y,n=oA n B n is the usual euclidean scalar product. The projection 



factor, which intervenes in the r.h.s. of (^) to keep b on a sphere : 

= <mr)Ar)> (23) 
< 6(r),6(r) > 

will be of central importance in the following. 

Characteristic time scales on shell n are in first approximation proportional to b n but now b n 



cannot exceed the initial value of y < b, b >. It follows that within the "projected dynamics" 
defined by Eq. ( p2[) , the cascade towards small scales is not accompanied by an acceleration of 
motion as in the original equations. There is now no impediment against taking a very large 
number of shells since the required time resolution does not grow anymore exponentially 
with N. By integrating numerically (p2|), we observed that any initial condition of finite 
support (i.e. b n (0) ^ for < n < n , with no N) gives birth at large times r to a 
solitary wave moving with a constant velocity towards small scales. In other words, a period 
T may be defined such that asymptotically, for r — ► +oo (a more precise condition reads 
l<r< NT, because some reflection will ultimately occur on the ultraviolet boundary), 

b n+1 (r + T) = b n (r) = b n (r) = b(n-^) (24) 
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Note that (^4]) implies A(t + T) = A{r). The shape of the final solitary wave is found to be 
always the same, up to the scaling symmetry 

b(n - 1) - Xb(n - A~) (25) 

and it is remarkably stable, as demonstrated by Figs. |9| and |K]. 

Let us now make the connection between this finding and self-similar solutions in shell 
models. This is easily done by writing any solution b{r) of Eq. ( p^ ) in the form 

6(r) = exp(- f A(r')dT') c(r) = B{t)c{t) (26) 

Since the non-linear kernel is quadratic, one gets for c(r) 

b(t)N[c] 

dr 

d _ 

The original dynamics : ~~^ = -^P] is recovered, after defining the physical time t as 

t( r ) = f T B(r')dr' (27) 

JO 

These straightforward manipulations prove that every solution b{r) of Eq. (|22|) can be 
mapped onto a solution c(t) of the real physical problem in the inertial range, according to 
the transformation law : 

c{t) = exp( / A(T')dr') b(r(t)) (28) 
J o 

where r(t) is obtained from the inversion of Eq. (p7|). It can now be seen that a travelling 
wave in the projected dynamics, of period T and average value in time (A(t)) > 0, is the 
signature of a self-similar solution in the true dynamics. Indeed, according to (|28|), each time 
the component of b of maximal amplitude moves from one shell to the next, c is multiplied 
by exp((A)T). From a comparison with the initial Ansatz (p0|), one gets : 

Q z = exp((A)T) 

or 
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2 - r£ < 29 > 

This formula allows one to obtain accurate estimates for z, since both quantities (A) and T 
are easily measurable (and their product is left invariant as it should by the scaling symme- 
try (E5D). 



The method was first used to compute z for various values of /3 in the ON model. Re- 
sults are summarized in table |, where a comparison between z — | and the asymptotic slope 
7 of the a p - curve is also made. We find a reasonable agreement between these last two 
quantities, in view of the comparatively large errors in the estimate for 7. It is important 

to realize that the existence of self-similar solutions has nothing to do with the presence of 

2 

chaos. In the ON model, the exponent z decreases gently from 1 to — , as p varies between 



0.145.. and 1 (as already noticed in |10|,|TT|, one has z = 1 for (3 < 0.145..). The analytic 
stucture of the solitary wave remains the same. The scaling function / of Eq. (f21~|) presents 
an essential singularity /(£) ~ 2^exp(— 2^) for £ — > +00 (up to subdominant multiplicative 
corrections), and an exponential tail /(£) ~ Q z ^ for £ — > —00. In the case where z = 1, the 
exponential tail is replaced by a second essential singularity /(£) ~ 2~^ exp(— 2""^). Table | 
shows that z ~ 0.88 at the transition between the regular and chaotic regimes, located near 
(3 = 0.349. This high value explains why the ON model (at least for Q = 2) is bound to 
exhibit rather strong intermittency in the chaotic part of its phase diagram. 
We were curious to extend this analysis to the complex GOY model. It is a simple matter 
to generalize Eq. ( ^2[) to the case of a complex vector. Details will not be given here. The 
conclusion of our (partial) investigations is that the GOY model also possesses only one 
ideal self-similar solution for a given value of e. Furthermore, this self-similar solution is 
purely real and positive, up to the phase symmetry (|12"D of the model. This means that the 
complex amplitudes 6„(r) take the asymptotic form b n (r) = e l9 ^b(n — ^), where the phase 
9(n), subject to the constraint 6{n) + 0{n + 1) + 0{n + 2) = 0, is the only footprint of the 
initial condition and the amplitude b presents a shape quite similar to the one obtained for 
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the ON model. Quantitative results are presented in the table |T[ The comparison between 
large order statistics and scaling properties of self-similar solutions was done only for two 
values of e : e = 0.75 and e = 0.5. While in the former case the same agreement is ob- 
tained as for the ON model, we find in the latter a discrepancy by a factor 2 between 7 and 
z — I . The discrepancy is even bigger if one extrapolates from the She-Leveque formula ( IB ) 
7 = 2/9 = 0.222. We think that the failure of self-similar solutions to explain intermittency 
at high orders in this case, lies in the closeness to the Kolmogorov value | of their scaling 
exponent z. After all, with small as 0.052 and a Reynolds number R = 10 as 

in our computations, the amplitude of singular fluctuations grows, upon propagating from 
the integral scale to the dissipative one, by a factor Q^ z ~^ nd = which does not 

exceed 1.5 ! This gives very little chance to such a fluctuation to survive collisions with 
the turbulent background and to govern statistics at large orders. We find it plausible that 
the mildness of singular fluctuations and the finite length of the cascade combine to pro- 
duce a new kind of intermittency with a more pronounced multifractal character. It is an 
interesting issue, left for further investigation, to understand how the system is then able to 
develop an asymptotic growth of the a p with p steeper than the one expected on the basis 
of self-similar solutions. 



V. MOVABLE SINGULARITIES AS A SIGNATURE OF CHAOS 

From a formal point of view, self-similar solutions studied in the previous Section de- 
scribe the approach of the system towards blowing-up, which, in the absence of dissipation, 
happens in finite time. It is also of interest in the context of nonlinear o.d.e.'s to consider 
movable singularities taking place at complex times. The local structure of such objects is 
intimately linked to the non-linearity, while their distribution in the complex t-plane may 
help to understand such physical properties as high-frequency intermittency Besides, 



according to Painleve's criterion, non algebraic singularities indicate usually lack of inte- 
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grability. This yields a very economical way to detect analytically the presence of chaos in 



any dynamical system (see for instance []T3J] and the references therein related to this topic). 
There are two main reasons why we report in this Section a study of movable singularities 
in the ON model, which at first sight is disconnected from the rest of the paper. The first 
one is purely technical : it turns out that the method used to determine the local structure 
of movable singularities (and possibly their position) is quite close to the one developed in 
Section 4 for tracking self-similar solutions. The second reason has more to do with physics : 
whereas self-similar solutions by themselves had nothing to tell us about the chaotic prop- 
erties of the model, we shall see that movable singularities in the complex t-plane disappear 
(or better said, get trapped on the last shells near the ultraviolet boundary) as (3 exceeds 
a value of order 0.394 ± 10~ 3 . It is tempting to speculate that this threshold marks the 
ultimate boundary between chaotic and regular dynamics, the one reached in the limit of 
infinite Reynolds number. We shall also get strong indications that movable singularities in 
shell models parameterize energy backflows and as such could be responsible for the peeling 
off of coherent structures as they cascade downwards to small scales. It will become rapidly 
clear to the reader that the analysis to be presented below, though restricted to the ON 
model, can easily be applied to any shell model with presumably similar conclusions at the 
end. 

We shall work with the vector b defined in Section 4. The quadratic degree of non- 
linearities implies that the only movable singularities are poles so that : 

b n(t) ~ for t -> t* and < n < N (30) 

where t* is an arbitrary complex critical time. The N + 1 residues a n form a vector a, which 
after substituting (|30| ) into ( f22"|) is seen to obey the condition : 

-d=N[d] (31) 

The problem now is to solve fl3~T|). This is a much more difficult task than computing fixed 



points as in Section 2. First, a is necessarily complex (it is easy to check that (pl|) implies 
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J2n=o a nQ~ 2n = 0)- Second, we expect on physical grounds the vector a to be localized in 
shell space. This means that we are looking for solutions of Eq. (|3l|) which would be square- 
summable (X}i£o | fl n| 2 < +°°)> were the range of shell numbers extended to the whole set 
of relative integers. Any "shooting" method of the type outlined in Section 2, which would 
start from one endpoint and try to join the other one with the appropriate asymptotic 
behaviour, is in fact doomed to failure because of strong numerical instabilities. 
As in the preceding Section the idea will be to approach dynamically the desired solutions to 



Before doing so, we must say a few words about the notion of "genericity" of movable 
singularities. Consider a singularity at time t* and assume a is known. Equation (|30f) gives 
only the leading order term in the expansion of b near t*, which may be pursued order by 
order just from local analysis. Writing b as - — — + 5b, where the correction 5b is small 
compared to the zero* h order term, one gets to linear order in 5b 

±5b = ^-M5b + F-°® (32) 
dt t-t* t-t* v ' 

In the equation above, M = -^-r— is the Jacobian matrix of the nonlinear kernel N\b) 

I dh 3 J 

evaluated at point a. Forcing and dissipation were kept for completeness in the right hand 
side of ( [52"D but only the homogeneous part of the equation really matters in what follows. 
It has (N+l) independent solutions of the form (t — where is the i th eigenvalue 

of M and 6j the corresponding eigenvector. Provided Re \Xi > — 1, a correction of the type 
\i(t — t*) flt bi with Aj an arbitrary complex number, is free to appear in the expansion of 
b around t*, since it is indeed smaller than the zero th term. Actually, N being quadratic, 
N[a] = —a implies Ma = —2a. Therefore, one of the (say /x ) equals by construction 
-2. The eigenvalue -2 and the corresponding eigendirection a are associated to the arbitrary 
position of t* and as such must be excluded from the expansion of 6. It follows that the 
most general expression of b around a singularity reads : 

a Ns 

b(t) = —— + J2 Ht - t*Y% + h.o.t. (33) 

t t i=l 

where N s is the number of eigenvalues of M, whose real part is bigger than -1. It is not 
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difficult to check that, once the N s complex numbers Aj are given, there is no arbitrariness 
left in the rest of the expansion (denoted as h.o.t. in fl33|)). What we have in our hands 
is a local expression of our solution which depends on (N s + 1) parameters (t*, X\, . . . Xn s ), 
whereas (N + 1) initial conditions are necessary to specify entirely the evolution of the dy- 
mamical system. Therefore a singularity will be generic (i.e. it will not result from a set of 
initial conditions of zero measure), if and only if N s = N. In other words, we are interested 
only in solutions to fl3lf ) with N eigenvalues \ii of real part bigger than -1, besides the trivial 
one /io = — 2. 



The previous considerations suggest the introduction of the following dynamics : 
d _ „-» r - .„ (a,N\a\) . rT (d,N\a])._ 

—a = -N[a\ + (Re x ' i J/ +iSIm x ' i u )a 34 
dr (a, a) (a, a) 

where, since we are dealing now with complex-valued vectors, (A, B) = Y^n=o ^n^n The 
second term in the r.h.s. of ( |3~4"D keeps the norm of a constant. The last one affects only its 
phase and one is in principle free to choose any value for the parameter 8. It may be shown 
that there is a one-to-one correspondance between fixed points of the dynamics ([53]) with 
a basin of attraction of finite measure, and generic solutions (in the sense of the previous 
paragraph) to the initial problem. A proof of this almost intuitive statement is given in 
the Appendix. It has nice consequences : in order to determine the possible arrangements 
of residues a n , it suffices to integrate (0) for initial conditions which are not purely real 
(otherwise they remain so forever). If after a long enough time, a stationary state df is 
reached, then : 

d=-d f (35) 
f (N[a f ],a f ) 

contains the desired information. Note that the computational cost of the method increases 
only linearly with the number of shells N + 1. It is therefore easy to get rid of finite size 
effects if necessary. 

We have applied this technique to the ON model and made the following observations. As 
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anticipated on the basis of the preceding considerations, the vector a evolves systematically 
towards a fixed point provided the condition 5 > 1 is met (actually, the marginal case 5=1 
still works but requires longer times of integration). After performing the rescaling (pop, the 
final state of a (giving access to the residues a n ) was found to be always the same, up to 
complex conjugation (which is an obvious symmetry of Q3T]) ) and translation along the shell 
number axis. This last property, which is crucial to ensure the "mobility" of the singularity 
in momentum space, holds for (3 < (3* = 0.394 ± 10 -3 . For (3 > (3* we find only one solution, 
rigidly attached to the last shell. Figures |Tl| and |12] summarize the phenomenon by showing 
the modulus and the real part of a n for respectively f3 = 0.39 and f3 = 0.40. They were 
deduced from a numerical integration of (|34]) with N = 29, 5 = 2 and the initial condition 
a n = i5 n $. The imaginary part of a n has not been represented in order not to burden the 
figures. For (3 = 0.39, a change in the initial conditions or in the value of 5 most likely leads 
to a displacement of the peak of the final structure along the horizontal axis. In contrast, for 
(3 = 0.40, the peak resides always on the last shell. A perfect convergence onto a true fixed 
point of fl3"4] ) is difficult to achieve because of slow transients near the transition. Thus we 
cannot exclude some minor adjustments of residues with respect to the picture shown here, 
especially at the rear end of the structure (n > 10 in Fig. |IT]). Note the characteristic pat- 



tern at the front (5 < n < 7 in Fig. |TT]) with a large negative excursion of Re(a n ), which by 



the way may be still recognized in Fig. [12], i. e. beyond the threshold. From a mathematical 
point of view, solutions in the inertial range, as depicted by Fig. [11], disappear when one of 
the eigenvalues /ij of the Jacobian matrix M gets a real part smaller than -1. Apparently, 
the only place where they manage to survive is near the ultraviolet boundary, where the 
nonlinear kernel is strongly modified. 

We shall not expand too much on these findings. At least they prove that complex time sin- 
gularities are not involved in the building-up of self-similar solutions, because in contrast to 
the former, the latter were found to exist for any value of f3. Just from this obvious remark, 
it is tempting to infer that complex time singularities in shell models, when sufficiently close 
to the real time axis, encode the occurence of "blockades" in the energy cascade, leading 
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possibly to negative excursions of shell amplitudes and more or less developed energy back- 
flows. This interpretation is corroborated by the wild oscillations displayed by the phase of 
residues and also the fact that such objects form most naturally at the ultraviolet boundary 



as suggested by Fig. 12. The system is bound to exhibit regular dynamics for (3 > (3* because 



it has lost these agents of disorder. 

Before closing this Section, we would like to mention that equations ([34]) , which were 
introduced as an abstract auxiliary tool, may also be used more concretely for locating sin- 
gularities of a real solution of the shell model (neglecting forcing and dissipation). Consider 
indeed initial conditions of the form a = ibo where bo is arbitrary but real. It may be checked, 
by using manipulations similar to those leading in Section 4 to Eq. (|28|), that stepping for- 
ward the fictitious dynamics ( p4[) is in fact equivalent to integrating the original dynamics 
—b = N[b] (from the initial condition bo) along a trajectory in time space parameterized 
as : 

t = -i / T exp( / T A{T")dr")dr' (36) 

JO JO 

where A(t) reads : 

A{T) = Be S^M. + a < 37 > 

{a(T),a(r)} {a{T),a(T)} 
For 5 = 0, A(t) is real and according to (|36|) the path followed in the complex time plane is 
parallel to the imaginary axis. The probability of crossing a singularity in this way is obvi- 
ously null for arbitrary initial conditions. For finite values of 8, the trajectory gets curved 
in such a way that, for 5 > 1, it finds with probability one a singularity of b(t) at the end. 



VI. CONCLUSION 



Starting from a numerical investigation of the ON model, we were led to identify ele- 
mentary bricks in its dynamics, which must exist more generally in any scalar shell model. 
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Interestingly enough, they appear to have rather constrained structures. Naturally the con- 
struction of a statistical theory from these deterministic objects remains a hard task. But 
we think that a precise knowledge of their properties may help to formulate new questions. 
For instance the discrepancy found in the case of the GOY model for e = 0.5 between the 
asymptotic growth of scaling exponents of statistical moments and the strength of extreme 
fluctuations is a puzzling fact, which clearly deserves further investigation. Another isssue 
concerns the selection mechanism of the scaling exponent z of self-similar solutions whose 
present understanding is still poor. One must remember that the method developed in the 
paper is in essence dynamical. We cannot therefore exclude the existence of a larger mani- 
fold of solutions, out of which only the element with the smallest z would be systematically 
observed. Clearly more mathematically oriented work would be welcome to elucidate this 
technical point, which may be of some physical relevance. 
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FIGURES 

FIG. 1. The amplitudes of the first two shells as functions of the parameter a describing the 
dissipative range. 



FIG. 2. Plot of the quantity F{a) defined in the text for Q = 2, R = 10 5 and (3 = 0.348. 
FIG. 3. K41-like static solution for the same value of parameters as Fig. ^. 

FIG. 4. Scaling properties of £ n)P for various p and (3 = 0.33. 

FIG. 5. Scaling exponents o v versus p for (3 = 0.28, (3 = 0.33, and (3 = 0.343. For the lowest 
value of (3 and p > 3, the cr p 's fall on a single straight line. There is more curvature for the other 
two values of /3, indicating multifractality. Nevertheless a linear growth is recovered for p > 8. 

FIG. 6. A zoom on the evolution of the cross-over region with (3. A particularity of the ON 
model compared to the GOY model is that it must get very close to the transition before showing 
a clear curvature of the a p 's. As a consequence, a p is very close to zero for p < 3. 

FIG. 7. Comparison of the statistics of the ON and GOY models far from the boundary 
between regular and chaotic dynamics. 

FIG. 8. Comparison of the statistics of the ON and GOY models in a regime where multifrac- 
tality is well established. 

FIG. 9. Amplitudes 65 (r), 612 (t), 619 (r) within the "projected" dynamics for f3 =0.335 and the 
initial condition b n (0) = S n o. The emergence of a travelling wave is clearly demonstrated. 

FIG. 10. Time evolution of the quantity A(t) defined in Eq. ( p3|) of the text, for the same values 
of parameters as in Fig. |9[ The small oscillations of A around its mean value in the asymptotic 
regime are due to the discreteness of the lattice. 
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FIG. 11. Modulus and real part of residues for = 0.39, as obtained from an integration of 
the fictitious dynamics (|34| ) with 5 = 2 and the initial condition a n = i5 n $ (N=29). 



FIG. 12. Same as in Fig. 11 but for = 0.40. In this case the same pattern is always observed, 



independently of the initial condition. 
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TABLES 






z 


Z 3 


7 


error 


0.15 


0.996 


0.329 






0.28 


0.921 


0.245 


0.24 


±10~ 2 


0.33 


0.889 


0.223 


0.213 


±10~ 2 


0.343 


0.881 


0.214 


0.20 


±210~ 2 


0.348 


0.878 


0.212 






0.7 


0.721 


0.054 




regular dynamics 


0.8 


0.692 


0.025 







TABLE I. Exponents z of self-similar solutions in the ON model for various values of (5. The 



left columns present estimates obtained from Eq. ( |29| ) after a numerical integration of Eq. (22). 
The last digit is given with an ±1 accuracy. The last two columns present data extracted from 
statistical analysis. The comparison between columns 3 and 4 show that scaling properties of 
self-similar solutions account in a satisfying way for large order statistics in the chaotic part of the 
phase diagram, even close to the transition where z remains rather big. 



e z z — | 7 error 

0.398 0.684 0.018 regular dynamics 

0.5 0.719 0.052 0.12 ±310" 2 

0.75 0.888 0.222 0.23 ±10~ 2 

0.8 0.946 0.279 

TABLE II. Same quantities as in table | but for the GOY model. One observes that z takes 
rather small values everywhere in the chaotic part of the phase diagram (e > 0.398). The disagree- 
ment between columns 3 and 4 for e = 0.5 is too large to be imputable to a lack of statistics. 
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APPENDIX: 



In this Appendix we establish the equivalence between stable fixed points of the dynami- 
cal system (|34|) introduced in the Section 5 and generic movable singularities in shell models. 
First, we observe that static solutions of (|34D , if they exist, are such that N[df] = Xdf where 

(df, N[df]) 

the coefficient of proportionality A = — — zr4 — obeys : 

{a f , a f ) 

A = Re(X) + i5Im(X) (Al) 

If 5 7^ 1, A is bound to be real. The case A = corresponds to df being an inertial fixed 

point (N[df] = 0). We discard this possibility (whose occurence would be easily identified in 

practice) and concentrate on the more interesting case of a finite value of A. Then a = — ^- 

X 

verifies as it should N[a] — —a. 

Now we must ask about stability properties of df. By linearizing the system of differential 
equations fl3^D around their fixed point, one finds the following evolution of small perturba- 
tions Sa : 

^Sa = X(M +l)Sa + ... (A2) 
dt 

where the terms hidden behind the dots are all directed in the direction of df and have been 
omitted for simplicity. As in Eq. (|32"D of Section 5, M is the Jacobian matrix of first order 
derivatives of the non linear kernel N evaluated at point a defined above. It appears that 
the space "transverse" to df belongs as a whole to the stable manifold of 3f, if and only if 
all the eigenvalues /ij of M for 1 < i < N have a real part bigger (resp. smaller) than -1 
with A negative (resp. positive). The second possibility would lead to a divergence of the 
trace of M in the limit N — > +oo, which contradicts the assumption of a finite norm for a. 
We are thus left with A < and by the same token N eigenvalues Hi of real part bigger than 
-1, which is nothing but the criterion for genericity established in Section 5. 
Finally, let us consider perturbations along the direction of df or a. Since the dynamics ( |34D 
preserves the norm of a(r), they reduce to phase fluctuations which may be parameterized 
as a(r) = e ie ^df. The phase 6{r) is found to obey the equation of motion : 
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^-6 = A(5-l)sin0 (A3) 
dr 

Therefore complete stability of Sf requires 5 > 1 (since A < 0), as announced in the main 
text of Section 5. 
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